Scaling of the glassy dynamics of soft repulsive particles: a mode-coupling approach 
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We combine the hyper-netted chain approximation of liquid state theory with the mode-coupling 
theory of the glass transition to analyze the structure and dynamics of soft spheres interacting via 
harmonic repulsion. We determine the locus of the fluid-glass dynamic transition in a temperature - 
volume fraction phase diagram. The zero-temperature (hard sphere) glass transition influences the 
dynamics at finite temperatures in its vicinity. This directly implies a form of dynamic scaling for 
both the average relaxation time and dynamic susceptibilities quantifying dynamic heterogeneity. 
We discuss several qualitative disagreements between theory and existing simulations at equilibrium. 
Our theoretical results are, however, very similar to numerical results for the driven athermal dy- 
namics of repulsive spheres, suggesting that 'mean-field' mode-coupling approaches might be good 
starting points to describe these nonequilibrium dynamics. 
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I. INTRODUCTION 

An assembly of hard spherical particles undergoing 
Brownian motion, if it can avoid crystallization {e.g. 
due to polydispersity), at a sufficiently high volume frac- 
tion undergoes a glass transition to an amorphous solid 
state In experiments, structural relaxation of col- 
loidal hard sphere systems stops near volume fraction 
tp w 0.60 0. The phenomenology of this so-called col- 
loidal glass transition is strikingly reminiscent of the 
molecular glass transition observed upon decreasing the 
temperature in glass-forming liquids. When Brownian 
motion is negligible, a hard sphere system undergoes in- 
stead a jamming transition near ip sa 0.64; it acquires 
rigidity by building a mechanically stable network of con- 
tacts between particles Q . This is most readily observed 
in granular materials. 

While the glass and jamming transitions of hard sphere 
systems have been widely studied for a long time, the 
study of dense systems composed of soft repulsive par- 
ticles is, by comparison, in its infancy. Colloidal parti- 
cles with tunable softness are now routinely prepared in 



the laboratory [J. 
iblc grains abound 




and examples of compress- 
1 1 01 ] - The rheological, struc- 
tural, or dynamical properties of both types of systems 
are currently actively studied experimentally by several 
groups 0, U H 0, H S 03 ■ This justifies current theo- 
retical efforts to understand the behaviour of dense as- 
semblies of soft repulsive particles both at finite temper- 
atures [ill, Ell EH, EI| , relevant for colloidal systems, and 
in the zero-temperature limit relevant for granular mate- 
rials [II, El, E3 EI El- 

In this paper, we combine the hy per- netted chain 
approximation of liquid state theory [20j with mode- 
coupling theory [2l| to analyze the equilibrium structure 
and dynamics of dense systems of soft particles with finite 



range, harmonic repulsion [15| : 

V(r < a) = e(l -r/af, 



(1) 



where e determines the strength of the repulsion, a is 
the particle diameter, and r the distance between two 
particles. Particles separated by r > a do not interact, 
V(r > a) = 0. The control parameters for this system 
are therefore the volume fraction tp = ira 3 p/6, with p the 
number density, and the ratio of the temperature and e. 
In the following we use reduced units: we give lengths in 
units of a and temperature in units of e. 

The system of harmonic spheres was originally intro- 
duced in the context of the zero-temperature jamming 
transition [l5l ]. More recently, its behaviour was inves- 
tigated at finite temperature using molecular dynamics 
computer simulations [ill, El, EH . To our knowledge, this 
model system was not studied theoretically in the context 
of liquid state and mode-coupling theories. However, its 
physics should be similar to a number of similar mod- 
els such as Hertzian spheres [22| or the Gaussian core 
model [HH!]. 

In both thermal and athermal contexts various scal- 
ing relations were reported for the dynamics of harmonic 
spheres in the vicinity of the glass or jamming transi- 
tions occurring in the hard sphere limit [ll|, 0, El, E3| • 
This type of scaling is usually rationalized by postulat- 
ing that harmonic spheres in fact behave as an 'effective' 
hard sphere system with a 'renormalized' volume frac- 
tion [ljj. However, the scaling formulae are typically 
established using largely empirical procedures. Our pri- 
mary goal in this work is to derive this scaling behaviour 
using a liquid state theoretical approach in order to put 
it on a firmer basis, at least for systems in thermal equi- 
librium. 

The paper is organized as follows. We present the- 
oretical methods which we use to obtain the structure 
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and dynamics of harmonic spheres in Sec. |H] We then 
present results for the equilibrium dynamics in Sec. IIIII 
We describe the phase diagram in Sec. |lVl In Sec. |V] 
we analyze dynamic heterogeneity using three-point dy- 
namic susceptibilities. We conclude the paper in Sec.lVIl 

II. THEORETICAL APPROACH 

Our theoretical approach consists of two steps. First, 
we use the hyper-netted chain (HNC) approximation of 
liquid state theory to obtain the static structure of the 
harmonic sphere system. Next, we use the predicted 
static structure factor as input of the mode-coupling 
equations to obtain time correlation functions of har- 
monic spheres and we use these functions to determine 
the dynamical behaviour. In the next two subsections we 
briefly describe the two elements of our approach. 

A. HNC equations and static behaviour 
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FIG. 1: Convergence of the numerical solution of the HNC 
structure factor to its T = hard sphere limit for q near the 
first diffraction peak at tp = 0.522. For a given number 2" 
of discretized values of r used to solve Eq. ((2j), there exists 
a temperature above which the VT convergence in Eq. is 
obeyed. The vertical line shows the lowest value of T used in 
this work. 



The HNC approximation [20j results in a closed equa- 
tion for the pair correlation g(r) of a fluid. It reads: 



g(r) = exp[-/3V(r) + g{r) - 1 - c(r)], 



(2) 



where /3 = 1 /T and c(r) is the direct correlation function 
defined through the Ornstein-Zernike equation: 

g(r) - 1 = c(r) +p f dr'c(\r - r'\)[g(r') - 1]. (3) 



We solve Eq. ([2]) numerically using an iterative proce- 
dure. Since the direct correlation function is smoother 
than the pair correlation function, the HNC equation is 
more easily solved in terms of c(r). Given the solution for 
c(r) after i — 1 iterations, Cj_i(r), we obtain the solution 
at step i as follows: 



/ \ FT » , , OZ „ / \ FT" 

Ci_i(r) > Ci-i{q) * gi-i(q) 



HNC 



c(r) -> Ci(r) = ac(r) + (1 - a)ci_i(r), (4) 



where the first and third steps arc Fourier transforms 
performed using fast Fourier transforms, the second one 
uses the Ornstein-Zernike relation ([3]), the fourth one 
uses the HNC closure relation @, and the last step in- 
volves combining the new direct correlation function c(r) 
with Cj_i(r) using a mixing parameter a. Convergence is 
achieved when the difference between a(r) and Cj_i(r) in 
Eq. becomes smaller than some prescribed precision. 

Since we want to access very low temperatures where 
the pair potential in Eq. ([T]) becomes equivalent to a 
hard sphere potential, some attention must be paid to 
the discretization scheme we use. Additionally, the fluid 
develops medium-range structure when density increases 
and/or temperature decreases, so that the cut-off in real 
space must be large enough. We used a real space cut- 
off L = 32, and discretized the interval [0, L] using 2™ 



points, which is convenient for the fast Fourier transform 
algorithm. To properly represent the low temperature 
behaviour, discretization must be accurate enough that 
the factor 



Y(r,T) = exp[-/?(l - r) 2 ] 



(5) 



is correctly described even at very low T. We note 
that, for a given discretization, there necessarily exists 
a temperature, T n , below which Y(r,T) is not well- 
described. It is easy to see that this temperature scales 
as Tji 2 ^ . 

This behaviour is ilustrated in Fig. [T] which shows 
how the structure factor S(q, T) of harmonic spheres at 
(p = 0.522 reaches its T — > value for increasingly finer 
discretizations (data are shown for q near the first diffrac- 
tion peak). To get accurate results down to T = 10 -8 , as 
shown below, we need to use n > 18. The value n = 19 is 
used throughout this paper. Additionally, since we need 
structure factors for temperatures and densities in very 
narrow ranges near the dynamic singularities, we had to 
pay close attention to the accuracy of the numerical so- 
lution of the HNC equation in order to resolve very close 
state points. 

In Fig. Q] we also show that difference between the fi- 
nite temperature S(q,T) and its hard sphere (T = 0) 
limit scales as VT. This behaviour can be derived as fol- 
lows. Since g(r < 1, T = 0) = for hard spheres we can 
decompose the difference S(q,T) — S(q,0) as follows 

S(q, T) - S(q, 0) = 4np f drr 2 ^^ 5 (r, T) 

Jo q r 

f 00 sin v 
+ 4np drr 2 ^^[g(r,T) - g(r,0)]. 
Ji qr 

Using the fact that in the low temperature T — > limit 
Y(r,T) becomes a narrow peak located at r = 1, it is 
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possible to approximate the first term in this expression 
by 

2n 3 ^ py(r = l,T = 0)^Vf. (6) 
9 

where y(r,T) = g(r,T)/Y(r,T) is the cavity function. 
Since both g(r,T) - g(r,0) and S(q,T) - S(q,0) are of 
the same order, we conclude that 

S(q,T)-S(q,0)~VT, (7) 

and this scaling is satisfied, sec Fig. [1] The \/T tem- 
perature dependence of the difference S(q,T) — S(q,Q) 
will play a crucial role in the next section, Sec. IIIII We 
will use it to motivate the dynamic scaling in the vicinity 
of the hard-sphere transition, within the mode-coupling 
approximation. 

In a recent molecular dynamics investigation [ll| , the 
temperature dependence of the energy density, e(T), was 
used to relate soft spheres to hard particles. We can 
use a similar reasoning to predict the low temperature 
behaviour of e(T). The energy density can be expressed 
in terms of the pair correlation function 

e(T) = 2np [ drr 2 V(r)g(r,T). (8) 
Jo 

The convergence of e(T) to its T = value e(T = 0) = 
can be estimated by making explicit the Y(r, T) factor 
in g(r,T) and using the cavity function. In this way we 
obtain 

e(T)«^py(r = l,T = 0)T 3 / 2 . (9) 

Molecular dynamics investigation [Tl| reported a power 
law behaviour of the energy density, e(T) ~ T^ 1 , with 
exponent p crossing over from p = 3/2 at low volume 
fraction to p ~ 1.3 at larger volume fraction. This in- 
dicates that the low temperature scaling regime © is 
not accessible at large densities in molecular dynamics 
simulations. 



B. MCT analysis and dynamic behaviour 

The mode-coupling theory (MCT) [2l| was originally 
derived to describe the dynamics of Newtonian systems 
[2l| . An analogous theory was later derived for Brownian 
systems [26| . Here we briefly present the latter version 
of MCT. 

The starting point of the theory is an exact equa- 
tion for the time derivative of the intermediate scatter- 
ing function F(q;t) in terms of the so-called irreducible 
memory function, 

d t F(q;t) = -^F{q-t)-^dt'M^{q;t-t')d v F{q;t'). 

(10) 



Here Dq is the diffusion coefficient of an isolated Brown- 
ian particle. Irreducible memory function M lrr (q;t) can 
be expressed in terms of a time-dependent four-point den- 
sity correlation function evolving with the so-called irre- 
ducible dynamics. The main approximation of MCT con- 
sists in factorizing this four-point function. In this way 
this somewhat mysterious quantity is reduced to a prod- 
uct of two intermediate scattering functions. After an 
additional technical approximation (which was indepen- 
dently shown to be quite innocuous) one arrives with the 
following expression for the irreducible memory function 

+ q • (q qi)c(|q - qi|)] 2 F( qi ;t)F(\q - qj|; t). 

Equations (j 1 Oti 1 1|) allow us to evaluate the time depen- 
dence of the intermediate scattering function. The only 
input required is the static structure factor S(q). It en- 
ters into Eqs. (jlOlllj) (note that c(q) = (1 - l/S(q))/p, 
from Eq. (|3|)) and it also provides the initial condition for 
the intermediate scattering function, F(q;t = 0) = S(q). 
It is easy to see that the natural time unit for our system 
of harmonic spheres is a 2 /D . In the following all times 
are given in terms of this unit. 

Numerical solution of Eqs. (|10H11[) is somewhat com- 
plicated because one needs to describe evolution of the 
intermediate scattering function on very widely separated 
time scales (see Fig. [2j . The commonly used algorithm 
was first described in Ref. [I?]]; here we use the imple- 
mentation described in considerable detail in Ref. [281 ] - 
Briefly, the basic steps to the algorithm are as follows. 
The intcgro-diffcrcntial equation is discretized and solved 
for 2iV s steps with a finite time step of St using any suit- 
able numerical algorithm. After 2N S steps are complete, 
the time step is doubled and the results from the initial 
2N S steps are mapped into a new equally spaced set of N s 
values for the quantities needed to continue the numer- 
ical algorithm. This mapping includes the integrals as 
well as the intermediate scattering functions. Then the 
numerical algorithm is restarted with the new time step 
and continued for another N s time steps and the map- 
ping is performed again. This procedure is continued 
until a convergence condition is satisfied. In the present 
work we used 300 equally spaced wave- vectors with spac- 
ing 6 sa 1.96, the first wave- vector at fco = 5/2 and the 
largest wavevector at A: max f» 58.81. 

In Fig. [5] we show the prediction of the mode-coupling 
theory for the time dependence of the normalized in- 
termediate scattering function F(q; t)/S(q) for harmonic 
spheres at ip = 0.5235. With decreasing temperature 
the relaxation becomes progressively slower. One should 
notice that an intermediate time plateau is developing 
which is the manifestation of the cage effect: with de- 
creasing temperature particles are trapped longer and 
longer within their first solvation shells and the final 
(a) relaxation shifts to longer and longer times. MCT 
predicts that at this volume fraction at temperature 
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FIG. 2: Time dependence of the normalized intermediate 
scattering function F(q;t)/S(q) as predicted by the mode- 
coupling theory for harmonic spheres at tp = 0.5235, for 
wavevector corresponding the first peak in the static struc- 
ture factor, q = q max ~ 7.17. The static structure factor 
used as input in mode-coupling equations was obtained from 
the hyper-netted chain approximation. The lines correspond 
to the following temperatures (from left to right): 0.0015, 
1.61 x 10 -5 , 5.39 x 10~ 7 , 1.00 x 10~ 7 , and 6.84 x 10" 8 . 



T c (ip = 0.5235) « 6.594 x 1(T 8 the plateau extends to 
infinite times and the system undergoes an ergodicity 
breaking transition which is commonly referred to as the 
glass transition. 

In this work we are primarily focused on the tempera- 
ture and volume fraction dependence of the a relaxation 
time which we define in the standard way, -Fornax! T a ) = 
e _1 , where q m&x _ is the position of the first peak in the 
static structure factor. MCT makes a number of detailed 
predictions for various aspects of the time dependence of 
intermediate scattering functions (including power law 
approach to and departure from the intermediate time 
plateau, and the wavevector dependence of the relaxation 
time) . The investigation of the temperature and volume 
fraction dependence of these properties are left for future 
work. 

The hard sphere glass transition within MCT is usu- 
ally investigated using as the input the static structure 
factor obtained from the Percus-Yevick (PY) approxi- 
mation. For the discretization used in this work MCT 
combined with the PY structure factor predicts the glass 
transition at (p^ ps 0.5159. Moreover, MCT predicts 
that upon approaching ipF c Y the a, relaxation time di- 
verges algebraically, r a ~ (Vc Y — V) r > with the 
exponent <-y HS > PY ps 2.59. 

In this work we use the hyper-netted chain approx- 
imation for the static structure factor because we are 
mostly interested in effects of finite temperature. How- 
ever, we anticipate that the low temperature results will 
be influenced by the behavior of the hard sphere sys- 
tem. Therefore, we also solved MCT equations in the 
T — > limit using as input the structure factor predicted 
by the hyper-netted chain approximation in this limit. 
For the discretization used in this work MCT combined 



with the HNC structure factor predicts the glass transi- 
tion at ip c ps 0.52315. Moreover, according to MCT upon 
approaching this critical volume fraction the a relaxation 
time diverges algebraically, 



with the exponent 7 ps 3.26. 

It is well known that in real colloidal systems the glass 
transition predicted by mode-coupling theory is avoided. 
Typically, as shown in recent contributions 0, EI HI, 
one can find an intermediate range of volume fractions in 
which the volume fraction dependence of the a relaxation 
time can be fitted to a power law with the exponent close 
to that predicted by MCT (one should note that when 
this procedure is used the critical volume fraction is one 
of the fitting parameters; typical values obtained from 
fits are about 10% different from MCT predictions). The 
result is that MCT-predicted power law divergence of the 
a relaxation time describes well the experimental and 
simulational data over approximately 3 decades of r a . 
Recent experimental and simulational results reported 
in Refs. 0, EH HfJ suggest that upon increasing volume 
fraction further this approximate power law is followed by 
an 'activated' regime according to which the a relaxation 
time has an essential singularity divergence at a higher 
volume fraction. 



III. DYNAMIC SCALING AT EQUILIBRIUM 

The tools described in the previous section allow us to 
obtain, for any given state point (ip, T) , the relaxation 
time T a ((/?, T) of the harmonic sphere system within the 
MCT approximation. We will report results for a broad 
range of volume fractions, ip <G [0.51,0.90], and tempera- 
tures, re [io- 8 ,io~ 2 ]. 

In Fig. [3] we show the evolution of r Q (ip, T) in the vicin- 
ity of the hard sphere glass transition occurring at T = 
and tp c ps 0.52315. The behaviour observed at finite tem- 
perature is easily explained. When ip < (p c , the relax- 
ation time increases when T decreases, but it saturates 
at low temperature to its hard sphere value which is finite 
at these densities. When increasing the volume fraction 
closer to ip c , this hard sphere value becomes larger, and 
the low temperature limit is reached at a lower temper- 
ature. For ip > ip c the system is a hard sphere glass in 
the T — ► limit and so T a (ip,T — ► 0) = 00. It is clear, 
however, that the system hits a finite temperature sin- 
gularity at a critical temperature, T c (tp) which increases 
continuously from T c (ip = tp c ) = when ip increases. 

At this stage of the description, these data resemble 
the ones found in numerical simulations [II], Q~l]. In 
Refs. [ll]> 03 a scaling analysis of the relaxation time 
was performed assuming that harmonic spheres at low 
temperature resemble an 'effective' fluid of hard spheres. 
Physically, this means that the softness of the potential 
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FIG. 3: Relaxation time of harmonic spheres as a function 
of temperature for various volume fractions in the vicinity of 
the hard sphere glass transition. The lines are the analyti- 
cal formula (|14[) using the scaling functions in Eq. (|20|l ; full 
lines correspond to (p < ip c , dashed lines to (p > tp c . The 
power law in Eq. (|17|l for ip = ip c w 0.52315 is shown with 
a dotted line. Corrections to scaling are seen for ip < 0.515 
and ip > 0.53. Volume fractions are: (i) 0.51, 0.515, 0.518, 
0.52, 0.521, 0.522, 0.5225, 0.5228, 0.523, 0.5231 (ii) 0.5232, 
0.52334, 0.5235, 0.5237, 0.524, 0.5245, 0.525, 0.526, 0.5275, 
0.53, 0.535. 



allows small overlaps between particles at low tempera- 
tures, so that the effective radius of the particles is re- 
duced by thermal fluctuations. While the temperature 
dependence of the energy density was used to estimate 
the average overl ap, and in turn, the effective hard sphere 
diameter in Ref. [Il[ , in the context of the mode-coupling 
approach a different route should be used to map soft to 
hard particles. 

Within MCT, at a fixed number density the dynamics 
of the system is uniquely controlled by the evolution of 
the static structure factor. This suggests that the most 
efficient way to map soft to hard particles is by matching 
the structure factor of harmonic spheres at {ip, T) to the 
one of a hard sphere system at an effective volume frac- 
tion tp e f[ = ip e g{ip,T). Combining the low temperature 
behaviour of S{q, T) shown in Eq. ((7|) with the fact that 
S(q) is smooth function of the volume fraction we easily 
find that 



ifieff{<P, T) S3 ip - CVT, 



(13) 



where c is a positive prefactor with subleading dependen- 
cies on temperature and volume fraction. 

The discussion in the two previous paragraphs leads 
to the following form of dynamic scaling that should be 
obeyed by the relaxation time of the harmonic sphere 
system: 



r a {p,T) 



TO 



\<Pc-V\~t 



hs fi 



\<Pc - <P\ 



T 



(14) 



where the scaling functions f±{x) respectively refer to 
volume fractions above and below tp c . To be consistent 




10 1 

T/\6<p\ 

FIG. 4: Collapse of the data shown in Fig. |3]using the scaled 
variables suggested by Eq. (| 14f) for all data in the range 
0.515 < ip < 0.53. This data collapse involves no free pa- 
rameter. The lines through the points represent the empirical 
scaling functions f±{x) in Eq. pO|l . 



with the qualitative behaviour described above, the scal- 
ing functions f±{x) must have the following limiting be- 
haviors: to recover the hard sphere plateau at low T 
below ip c we need to have 



f—{x — > oo ) ~ const; 



(15) 



whereas to have a well-behaved r a {ip,T > 0) across ip c 
we have to require 



f-{x -> 0) ~ f+{x -> 0) ~ x-< 
The latter limiting behavior implies that 



T a {ip = <p c ,T)~ [ — 



7 HS /2 



(16) 



(17) 



Finally, a finite temperature algebraic singularity is ob- 
tained for ip > p c if there exists some x± such that 



f+{x -> X*) ~ (.T* - x)~ 



(18) 



We have found excellent agreement of the MCT pre- 
dictions with the scaling form in Eq. (|14[) . which sug- 
gests that the relaxation times at various ip and T can all 
be collapsed along two branches by plotting the rescalcd 
time, |<5y| 7 r a {ip,T), as a function of the rescaled dis- 
tance to the critical point, \Sp\/VT, where 8ip = ip c — p. 
The data collapse is presented in Fig.|U It works remark- 
ably well for the range of volume fraction 0.515 < p < 
0.53. It should be noted that no free parameter is in- 
volved in this data collapse, which is uniquely controlled 
by the T = hard sphere results from Eq. (|12jl . and 
by using the appropriate relation between temperature 
and density from Eq. (| 13|) . In Refs. [ll|, dS]' a similar 
data collapse was used in the opposite direction to infer 
the hard sphere behaviour from the finite temperature 
dynamics of the harmonic sphere system, a philosophy 
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which is clearly supported by the results presented in 
this section. 

The next step is to take Eq. (fl3|) more literally and to 
combine it with the results for the dynamics of the hard 
sphere system discussed in Sec. lIIBI to make the following 
ansatz 

T a (tp,T) K 7% S [tp eS (tp,T)}. (19) 

Eq. 1X9]) leads to the following scaling functions: 

/iW = (-T«)" /S , (20) 

x 

with a « 0.72 being the only adjustable numerical fac- 
tor (note that a is related to c in Eq. (fT3")) ) since the 
values 7 HS = 3.26 and tp c — 0.52315 are directly taken 
from hard sphere results, while the scaling variable x = 
\Stp\/VT was derived in Sec. Ill Al Clearly, these scaling 
functions are fully compatible with the constraints de- 
scribed in Eqs. |T5l QH HHJ)- The scaling functions (|20]) 
are shown as lines in Figs. [3] and 2] 

We have noted several times the similarity between 
the present theoretical results and the numerical results 
and analysis in Refs. [ill, E3- I n the simulations, a scal- 
ing analogous to the one in Fig. |4] was presented for the 
behaviour of \ogT a (tp, T) instead of T a (tp,T) here. This 
implies that the scaling behaviour predicted by MCT is 
in fact in strong quantitative disagreement with numer- 
ical results. This should not come as a surprise since 
MCT is not able to describe the thermally activated re- 
laxation which takes place in real glass- formers. Thus, 
MCT predicts algebraic divergences which are never ob- 
served in simulations and experiments, and are replaced 
by stronger, generically exponential, divergences. 

It is interesting to note that algebraic scaling be- 
haviours and divergences seem to be well obeyed in the 
case of non-equilibrium driven athermal dynamics stud- 
ied in Refs. [TBI, [H, Oil ■ This is again not surprising since 
in these dynamics the system simply relaxes to the near- 
est energy minimum without being able to cross energy 
barrier using thermal activation [30| . This suggests that 
mode-coup ling approaches and 'mean-field' models (see, 
e.g. Ref. [3l|) might well be excellent starting points 
to tackle the athermal driven dynamics of soft repulsive 
spheres. 

IV. PHASE DIAGRAM 

In this section, we move from scaling properties very 
near tp c and give a broader perspective on the behaviour 
of the system in the (tp, T) phase diagram. The scaling 
results presented above suggest that the system is ergodic 
at all temperatures when tp < tp c . For tp > tp c , Eq. (|20p 
predicts that the relaxation time diverges at t q ~ (T — 
T c (tp))~' y with a critical temperature which vanishes 
continuously at tp c as: 

T c (tp>tp c )~(tp~tp c ) 2 . (21) 
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FIG. 5: Phase diagram obtained from the MCT analysis. 
Filled circles are transition temperatures obtained in this 
work, which follow the scaling behaviour in Eq. pip in the 
vicinity of tp c , as shown by the full line. Open symbols are 
the mode-coupling critical temperatures obtained in Ref. 

Since the scaling behaviour in Eq. (TFf| only holds up to 
tp ~ 0.53 we have fitted the temperature dependence of 
the relaxation time at larger volume fraction to the usual 
power law divergence found within MCT, 

T a (<p,T)~{T-T D {<p))-*M, (22) 

using T c (tp) and "f(tp) as fitting parameters. 

In Fig. [5] we show the evolution of the resulting T c (tp) 
which thus delimits the fluid and glass phases in the the- 
oretical phase diagram of the system. These data con- 
firm that the scaling behaviour in Eq. (|2~lj) of the criti- 
cal temperature is only obeyed in the vicinity of tp Cl and 
clear deviations are seen at larger volume fractions where 
the scaling prediction overestimates T c by quite a large 
amount. 

We also find that increasing the volume fraction affects 
the value of the critical exponent 7. While 7(93 w tp c ) = 
7 HS w 3.26, we find that 7 decreases rapidly when tp in- 
creases: 7(0.58) « 2.92, 7(0.62) w 2.71, 7(0.68) w 2.50, 
7(0.75) w 2.41, and 7(0.90) « 2.37 This is a clear indi- 
cation that when moving away from tp c the system also 
leaves the universality class of the hard sphere transi- 
tion. This means that it becomes impossible to describe 
the soft spheres as 'rcnormalizcd' hard spheres when tp 
becomes too large. 

We confirm this statement in Fig. [5] where we show the 
evolution of the static structure factor along the critical 
line T c (tp). While nearly perfect collapse of the data is 
obtained for tp c < tp < 0.53, small deviations become 
noticeable for tp « 0.535, and are considerably amplified 
when tp increases further. It is this large difference in 
the shape of the structure factor at the critical temper- 
ature which accounts, within MCT, for the continuous 
evolution of the critical exponent 7. 

We can again compare these theoretical predictions to 
the numerical analysis reported in Ref. [l 21 ] . Although 
the mode-coupling algebraic singularity is not observed 
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tp c < ip < 0.53 
0.535 <ip< 0.75 




FIG. 6: Evolution of the static structure factor along the 
MCT critical line T c (p>). Full lines show that volume frac- 
tions between ip c = 0.52315 and ip = 0.53 collapse on the 
hard sphere structure factor at ip c , while deviations appear 
and increase rapidly at larger ip showing that soft spheres are 
not simply 'renormalized' hard spheres far above ip c . Volume 
fractions are: (i) 0.5234, 0.5235, 0.5237, 0.524, 0.5245, 0.525, 
0.526, 0.5275, 0.53. (ii) 0.535, 0.58, 0.62, 0.68, 0.75. 



in real liquids, it is usually found that such a power law 
behaviour is obeyed over a limited time window of ap- 
proximately 3 decades, which allows a rough determina- 
tion of the location of the 'avoided' mode-coupling sin- 
gularity. The outcome of this exercise for the dynamics 
of harmonic spheres as determined in computer simula- 
tions reported in Ref. [l2j is shown in Fig. [5] with open 
symbols. The mode-coupling line determined numeri- 
cally has qualitatively the same behaviour as the theo- 
retical line. Although the data do not span a very large 
temperature window, they are indeed compatible with a 
power law scaling as in Eq. (|2ip near the hard sphere 
mode-coupling singularity ip Cl and the critical line be- 
comes smaller than the scaling prediction at larger den- 
sity. However, the theory is quantitatively inaccurate as 
it significantly overestimates the critical temperatures at 
all (p. Note that for hard spheres MCT combined with 
HNC structure factor underestimates ip c by an amount 
comparable to that reported for MCT combined with PY 
structure factor. These discrepancies are well-known fea- 
tures of the mode-coupling approach [2l|, and MCT de- 
scriptions of real data usually imply analysis of scaling 
behaviour near singularities whose locations must be self- 
consistently determined by fitting. 

A more surprising disagreement between theory and 
simulations is the evolution of the critical exponent with 
density. While theory predicts a substantial decrease of 7 
at large volume fraction, numerical results indicate that 
7 increases instead very rapidly with ip above the hard 
sphere value [l2||. It is not clear whether this disagree- 
ment stems from an incorrect prediction of the liquid 
structure by the HNC closure, or from the mode-coupling 
approach itself. 



V. DYNAMIC HETEROGENEITY 

A newer and lesser known application of MCT is to 
use it to estimate the strength of dynamic heterogene- 
ity accompanying the glass transition using multi-point 
dynamic susceptibilities 32fl . It is well-known that dy- 
namics near the glass transition is spatially heteroge- 
neous, meaning that different parts of the system relax 
at different rates, while relaxation is correlated over a 
lengthscale which increases when the glass transition is 
approached [33| . 

A useful tool to quantify the strength of dynamic het- 
erogeneity is the four-point dynamic susceptibility Xi(t) 
which is defined from the spontaneous fluctuations of 
time correlation functions [H, KiH : 



Xi{t) = N[(f 2 (q-t))~(f{q;t)f 



(23) 



where f(q; t) represents the instantaneous value of the in- 
termediate scattering function F(q;t). Intuitively, Xi(~t) 
increases if correlations within the system get large, as 
the number of independently relaxing units within the 
sample decreases [3a]. Formally, Xi($) is a l so tri e volume 
integral of a spatial correlator quantifying the extent of 
correlations between local, spontaneous fluctuations of 
the dynamics and can thus directly be considered as a 
proxy for the number of particles that relax in a corre- 
lated manner close to the glass transition [37l |. 

We build on the results of Refs. [H, and esti- 

mate the spontaneous dynamical fluctuations quantified 
by X&{fy using linear response theory: 



, s T 2 fdF(q;t) 



(24) 

This relation is known to be an accurate representation 
of X4 (t) within the MCT approach [3^| and amounts to 
measuring the response of the averaged dynamics to ex- 
ternal fields in the linear regime (38j . 

The expression in Eq. (j2~4")) is highly convenient in the 
present context as we can directly obtain analytical re- 
sults for the scaling behaviour of X4(*) m the vicinity of 
(p c using results from the previous sections. Since we are 
interested in the scaling properties of the dynamic sus- 
ceptibility, we make two further approximations to obtain 
an analytical form. We first use the fact that the time 
decay of the intermediate scattering function obeys time 
temperature superposition, F(q;t) ~ !F(t/T a ). Thus we 
have (with x = T,cp): 



dF(g-t) _ t T , ft \ dUj a ) _ 
dx r a \T a J dx 



which is a non-monotonic function of time with a maxi- 
mum for t « T a . We focus on the height of this maximum, 
Xi = Xi(t — T a), which can then be estimated from the 
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FIG. 7: Evolution of the peak of the dynamic susceptibility 
X4 estimated from Eq. (|27[) for the same volume fractions of in 
Fig- EI including data below (full lines), at (dotted line), and 
above (dashed lines) <p c . The behaviour is clearly reminiscent 
of the one of r Q (ip, T) in Fig. [3] 



FIG. 8: Evolution of the contribution of the term contain- 
ing the temperature in Eq. (|24[) for the same parameters and 
using the same representation as in Fig. [7] Note that the ver- 
tical scale in both figures is different, and that the thermal 
contribution to X4 is always much smaller than the density 
contribution. 



behaviour of the relaxation time alone: 



cally: 



X4 = C 2 
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<91nT Q 
dT 



S(0, I> 2 



(91nT Q 

dip 



(26) 

with C = ^"'(1), so that C = (3e 1 for a stretched expo- 
nential lineshape, F{x) ~ exp(— x 13 ). 

To proceed analytically we make use of the scaling form 
in Eq. (Q3J) for the relaxation time. We get: 



X4 /C 2 = -^-G 2 ± (x)+S(0,T) 
4cy 



7 H V ! 



(1-G ± (x)/ 7 HS ) 2 , 
(27) 

where G±(x) = xf±{x)/ f±{x). 

In Fig. [7] we show the evolution of ^4, evaluated from 
Eq. ([2T)l across <p c . It is clear from this figure that \4 
has scaling properties very similar to the ones of the re- 
laxation time (Fig. [3]). It increases and saturates to a 
plateau when T decreases for ip < ip c , obeys a power 
law behaviour at ip c , and diverges algebraically at T c (<p) 
above ip c . 

In particular, we find that the term proportional to 
Cy 1 and stemming from the temperature derivative in 
Eq. (f2~4"|) is always safely negligible to the volume fraction 
derivative term in the dynamic range shown in Fig. [7] In 
fact, this figure would be almost unchanged if we had 
shown only the second term in Eq. (|2"4")l . Physically this 
implies, not too surprisingly, that dynamic heterogene- 
ity in the scaling regime of the harmonic sphere system 
is mainly controlled by density fluctuations, just as for 
hard spheres Q , while energy fluctuations play little role. 
We note that the opposite is true in supercooled liquids, 
where density fluctuations seem to be generically domi- 
nated by energy fluctuations [53, Hil . 

A second interesting consequence is that the scaling 
behaviour of x& near tp c can then be obtained analyti- 



VT 



where 



xf s (v) ~ <P 2 /(<Pc 



(28) 



(29) 



is the hard sphere result [3!|, and X±[x) = (1 — 
G±(a:)/7 HS ) 2 . The scaling behaviour in Eq. ([28]) is simi- 
lar to the one of T a (ip,T) found in Eq. (fT3)l . In fact the 
similarity is even quantitative, since combining Eqs. (|14l 
1201 [28)) , we can explicitly show that the relationship be- 
tween the four-point susceptibility and the averaged re- 
laxation time is identical for soft and hard spheres in 
the scaling regime near ip c , up to subleading contribu- 
tions. This suggests that plotting Xiifi T) vs. r Q (<^, T) 
would collapse all data for harmonic spheres onto the 
hard sphere data. 

For completeness we also show the 'thermal' contribu- 
tion to X4(t) m Eq. |24|) in Fig. because this term is 
too small to have observable effects in Fig. [7] Although 
its shape seems similar to the one of %4, it is not quite the 
same: it vanishes as T — > for ip < ip c , because T a does 
not depend on T in this limit. It follows a power law be- 
haviour for tp = ip c , but this divergence is in fact entirely 
due to the 1/cy prefactor, since the specific heat behaves 
as c v {T) - VT from Eq. ©. Finally, for ip > tp c both 
terms contributing to Xi diverge in the same power law 
manner, as (T — T c (ip))~ 2 , but the respective amplitude 
of the two terms is set by (ip/Sip) 2 and 1/cy- This im- 
plies that as long as ip is close to ip c the density derivative 
term dominates over the temperature derivative term. It 
is only when ip is much larger than tp c that the tempera- 
ture contribution might become dominant, but Xi is n °t 
described by Eq. (|27[) anymore and a direct evaluation 
of all contributions would be required to investigate this 
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crossover at large volume fractions. We have not pursued 
these investigations. 

To the best of our knowledge, dynamic heterogeneity 
has not been discussed numerically in harmonic spheres, 
and we cannot compare the present results with numer- 
ical results. However, for hard spheres, it has been es- 
tablished that the power law scaling in Eq. (|29|) is barely 
visible on actual data Q , and dynamic correlations seem 
to increase much more slowly with increasing the density 
than predicted by MCT, as is found also for supercooled 
liquids (37|. 

We note again, however, the close similarity between 
the present MCT results for dynamic correlations and the 
scaling properties found in driven athermal simulations of 
harmonic spheres where algebraic divergences of spatial 
correlations of particle dynamics and scaling properties 
very similar to Fig. [7] were reported pH fl9l. 1411] . 



VI. DISCUSSION 

In this paper, we have investigated theoretically the be- 
haviour of dense assemblies of harmonic spheres at low 
temperatures in a broad range of volume fractions, en- 
compassing the glass transition of hard spheres at ip c . 
We have combined hyper-netted chain closure for the 
structure with mode-coupling theory for the dynamics. 
We find that for finite temperatures near <p c , harmonic 
spheres behave effectively as hard spheres with a renor- 
malized volume fraction. This directly implies a scaling 
form for the relaxation time in the part of the volume 
fraction - temperature phase diagram in the vicinity of 
the hard sphere transition, which applies also to the am- 
plitude of dynamic heterogeneity. At larger volume frac- 
tion, deviations from hard sphere behaviour arise, and 
dynamic scaling breaks down. 

When compared to numerical simulations of the dy- 
namics of harmonic spheres at thermal equilibrium, the 
known shortcomings of mode-coupling predictions clearly 
show up: MCT predicts algebraic singularities that are 
not observed in simulations, and fails to predict the 
'activated' scaling behaviour observed in simulations of 
both hard and soft particles. In particular, using a 
mode-coupling approach we cannot discuss the large 
change of glass fragility with volume fraction discussed 
in Refs. [ll|, Q~l], as these come from subtle deviations 
from an Arrhenius behaviour which is not predicted with 
MCT. Another feature which is not well captured by the 
present calculations is the volume fraction dependence of 
the MCT critical exponent which increases with ip in the 
simulations, but decreases in our calculations. It is not 
clear, however, whether this last failure originates from 
the approximation used for the structure factor or from 
the mode-coupling theory itself. 

Although MCT is qualitatively unable to describe the 
nature of the glass transition in harmonic spheres, there 
seems to exist a time window of approximately 3 decades 
where its predictions can be applied. In experiments with 



hard sphere colloids, it is in fact only very recently that 
deviations from MCT behaviour were unambiguously ob- 
served in experiments covering a very broad range of re- 
laxation times [2] • This means in turn that for 'standard' 
experiments the behaviour predicted in the present arti- 
cle might still be of some value, and our theoretical ap- 
proach could certainly be extended to a broader family of 
soft pair potentials beyond harmonic interactions. Thus, 
we hope that the present results will motivate further 
analysis of the dynamics of soft colloids in simulations 
and experiments. 

Although we mentioned in the introduction that soft 
colloids are currently studied by several groups, the glass 
transition of soft colloids has only very recently been 
studied in a system made of microgel particles (7j. In 
this paper, three types of particles with increasing soft- 
ness were studied. The most striking result of this study 
is a change of the volume fraction dependence of the re- 
laxation time with softness, from r Q ~ (ip c — </2)~ 7 for 
hard particles, to log(r Q ) cx ip for very soft particles. 
Since the interparticlc interaction in this system is not 
known, one could imagine using a potential such as in 
Eq. |T]), with increasing temperature playing the role of 
increasing particle softness. Our results in fact predict 
that the volume fraction dependence of the relaxation 
time does not vary from the hard sphere behaviour for a 
very broad range of temperatures of at least 6 decades, 
see Fig. [31 Since the qualitative behaviour found in this 
work should be independent of the details of the pair 
potential, this suggests that, very likely, the change of 
particle softness in Ref. Q also corresponds to a change 
of the form of the interaction between particles, and is 
thus difficult to explain theoretically on the basis of the 
present work. 

Does the scaling behaviour found for harmonic spheres 
teach some lessons for understanding the glass transition 
of molecular glass-formers? A major conclusion drawn 
from the present theoretical results is that the physics 
of harmonic spheres is simply the one of hard spheres: 
they undergo a glass transition both upon compression 
or upon cooling with a relaxation time which diverges 
when the effective volume fraction tp e f{(tp, T) of the har- 
monic spheres system becomes equal to tp c , the hard 
sphere critical packing fraction. Moreover, the physics 
for <p ff(<p,T) < ip c is the same for both soft and hard 
particles, as seen for instance from the behaviour of the 
dynamic susceptibilities discussed in Sec. |Vj 

For this physical behaviour to be useful to understand 
aspects of the glass transition, one should invoke the pos- 
sibility that real liquids can be effectively described as 
hard spheres, an assumption which has a long history in 
the field of liquid state theory [2(| • It was revisited very 
recently in the present context in Ref. [13] which estab- 
lished that the dynamics of real liquids, and in particular 
the interplay between density and temperature, is qual- 
itatively different from the one of soft particles with a 
finite interaction range such as harmonic spheres. For 
instance the large change of glass fragility observed for 
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harmonic spheres is not observed in molecular liquids, 
which obey a much simpler scaling yielding glass fragili- 
ties independent of the density [42]. Here, we also found 
that the main contribution to the dynamic susceptibil- 
ity in Eq. (|24[) is always given by the density contribu- 
tion, while in real liquids the temperature term domi- 
nates [37], HO, H3| , suggesting that a different physics is 
at play in both cases. It could be, for instance, that at- 
tractive forces not included in potentials such as Eq. (p} 
provide a non-negligible contribution to the energy barri- 
ers that need to be crossed during structural relaxation. 

Finally, we comment on the intriguing similarity em- 
phasized throughout this paper between the present re- 
sults and the dynamic scaling behaviour discussed in sev- 
eral recent articles dealing with the athermal, driven dy- 
namics of harmonic spheres [HI, [10, |4l| . In both cases, al- 
gebraic divergences of dynamical quantities are obtained 
in the hard sphere limit, with a dynamic scaling be- 
haviour observed in its vicinity both at the level of the 
averaged dynamics and of the dynamical fluctuations. 
Note, however, the very different nature of the critical 
density in both cases [3(|: the hard sphere glass transi- 
tion discussed here is defined from the divergence of an 
equilibrium quantity, which is thus by definition indepen- 
dent of the preparation protocol of the system. Instead 
the zero temperature jamming transition does explicitly 
depend on which ensemble of configurations is selected 
by the studied dynamics [H, HH, [43|, There is no 

limit where these two distinct transitions can merge. 

We believe that the deep underlying explanation of this 
similarity is the fact that in both cases, the dynamics is 
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